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@ Methods and apparatus tar efHolent allocation of resources by optimizing nonlinear, convex functions with linear 
constnalnts* 

A method and apparatus Is descrltied for optimally 
focatlng resources. The optimal allocation Is done by 
minimizing a cost (which ia a convex non^near function of 
vBjIous aOocatlon variables) subject to cSfferent constraints 
(which are Gnear functions of the allocation variables). The 
method Mttally piclcs a state of the above variables (xo) in the 
Interior of the solution polytope (where the constraints are 
satisfied) and computes successive states xi^,« which 
progressiveV reduces the cost of allocation. The akwve 
iteration stops when suitable stopping ailes are met 

The method employs (i) an affine scaling transformation (a 
variant of t<armarkar's projective transformation) of the linear 
constraints, (ii) an ellipsoid to sphere transformation of the 
curved cost surfaces, (ill) a potential search scheme on the 
curved constant cost surfaces, (iv) an affine scale ad|ustment 
mechanism, and (v) a lir>o-eearch scheme. 
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Descripti n 



METHODS AND APPARATUS FOR BFFICIEMT ALLOCATION OF RESOURCES BY OPTIMIZINQ 
N0N4JNEAR. CONVEX FUNCTIONS WITH UNEAR CONSTRAINTS 



S Background of the Invention 

This Invention relates to systems for optimally aHocatfng resources among a plurality of resource utilizers. 
More particularly, this Invention relates to methods and apparatus which perform improved allocation of 
resources in a physical system whose allocations costs follow a non^inear function and whose allocated 
parameters are linearly constrained. . . ^ . , ^ 

10 The need for resource allocation and optimization of systems arises In a broad range of industnal and 
technological areas. Examples of such are assignment of transmission facilities In telephone transmission 
systems, control of product mix In a factory, allocation of people to do a specified set of jobs such as operator 
sendees provided by telephone companies. Inventory control, deployment of industrial equipment, and many 



others. 

IS • Resource allocation and optimization decfetons are typically subject to certain constraints on vanous 
decision variables. Resources, for example, are always limited in availability, and sometimes their usefulness in 
a particular application also may be limited. Thus, the problem is to allocate the resources such that the 
constraints are satisfied while at the same time the cost of allocation Is minimized, or the benefits arising from 
the allocation are maximized. The cost to be mmimlzed can be any measure selected by the user, such as 
20 minimizing expenditures, minimize delay, minimizing inefficiency, etc. . . ^ , 

A general method of characterizing and solving such optimization problems Is via mathematical 
progrBmming models. That Is, the various resources of the system to be optimized are represented by a 
variable (one variable for each resource assignment), the collection of variables is viewed as a vector, a 
solution for the desired state of the system to be optimized Is derived by evaluating a value for the vector 
25 through some selected procedure, and, lasUy. the values of the components of the evaluated vector are said to 
represent the desired values of the allocated resources. 

There are many physical optintization problems that, when modeled as described above, consist of linear 
relationships and thus belong to a class of problems and their solutions that are known as Unear Programming 
(LP) However there is also quite a large number of other physical optimization problems that, when modeled 
30 as described above, comprise some linear relationships and some non-liner relationships. These are not LP 
problems. For example, the total amount of cargo carried by an air canrter cannot exceed the capacity of the 
planes which the can-ier owns (at one cost) and which the carrier can lease (at different costs). This situation is 
typically modeled mathemailcaliy by a set of Inequality constrahts. and in such a situation, the cost function to 
be minimized may very well be non-linear. For example, the cost of leasing planes will go down (per ptene) as 
SS the carrier Increases the number of leased planes, until the supply of leasable planes diminishes to the point 
where per plane cost of leasing will begin to rise again. Such a relationship is called convex. 

In the description that follows, physical systems where the resource aUocatlon effort reajlts hi a 
mathematical model where the operating parameters are Hnearty constrained and ^^e wst ftmcton to be 
minimized is non-linear, are referred to as "Non Linear Functions - Unear Constrains (NLFLC) systems. Of 
40 particular importance is a subset of NLFLC systems which is characterized by a cost function that is quadratic. 
A quadratic cost function is a function that Is convex, with the additional property thai a line conrwcting any 
two points on the function do not cross any other point on the function. This is Bluslraledin RG. 1 • ^"^e 
10 is shown as a function of the two variables x and y. Ariy two points on curve 10. such as points 17 and IB. 
can be connected with a Dne 19. and line 19 does not Intersect cunw 10. Cuwe 10. then. Is quadratic. Curve 40. 
45 on the other hand, is not quadratic because points 39 and 38 are connected by line 46 and line 46 Intersects 
curve 40 at two other points. However, line A7 which connects point 39 with point 37 does not intersect curve 
40 at any other point. It can be said, therefore, that curve 40 has quadratic regions. x^^^,^«, 
Resource allocation tasks that can be modeled as NLFLC problems can be visualized in a manner WenttMl 
to that of conventional linear programming (LP) tasks as described, for example, seminal inwntiorj of 1^ 
50 Karmaricar. U.S. Patent 4744.028. issued May 10. 1988. To wit. each parameter which can be altered In the 
process of allocating resources is viewed as a dimension in a multi-dimensional orthogonal space. 
of feasible solutions Is found within the volume defined by the linear constraints of the system, each or wmcn 
describes a plane in the multi-dimensional space. The collection of the intersecting muttWimenslonal planes 
- forms that volume. The multl-dlmenslonal surface fomied by the intersecting muttl-dlmensional planes Is dteri 
55 called a -polytope". and the sides are often called •walls*, it has been known that, In general, the opwnaJ 
solution to a NLFLC problem lies on a wall of the pdytope. This is in contrast with the LP case (where the cost 
function Is also a linear hjnctton of the parametera) where the optimal soiutten fies on a vertex of the 

summary of prior art for solving NLFLC problems can be found in the book "Practical Optimization- 
$0 by P. E. Gill. W. Murray, and M. H. Wright, published by Acad mic Press In 1981. However, these metlwte 
simply cannot solve very larg problems, and the proW ms that these methods are capable of solying tal^ a 
relativety long time to reach a solution. In short none of these methods are suitable for solving even 
moderately-steed resource allocation prebl ms In real-time (that Is. sufficiently fast) to provide more or less 
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Brief DescripHon of the Drawing 
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continuous control of an ongoing process, system, or apparatus. The deficiency stems from the fact that these 
methods consume a tot of computer time before th y come up with an optimal allocation strategy. 

To overcome the computational cfifficulties associated with the known methods to solve linear programs. N. 
K. Kamiaricar frwented the new method disclosed in the eforementioned U.S. Patent No. 4.744.028. According 
to this method, a starting feasible point is selected from the Interior of the polytope. and a series of moves are 
made in th direction that, locally, points in the direction of greatest change toward the optimal vertex of the 
polytope. This direction, which is the direction of steepest •descent' does not depend on the parameters 
(variables) to be optimized. A step of computable size is then fallen in that direction, and the process is 
repeated until a point is reached which is close enough to the desired optimal vertex. One specific aspect of 
Karmari<ar's algorithm is that the above mentioned procedures are conducted in a transformed space and the 
outcome is mapped back into tha original space at each Iteration, it Is Importoit to remember that the 
transformation performed by ICarmarkar (called the "prolectlve transfonmatlon*| Is with respect to the 
constraints of the problem. No special transformation is done (as It is not necessary) with respect to the cost 
function. Under projective transfonmatton. the llneartty of the cost function is preserved. 

Alas, this method is not suitable for optimizing general Nli=LC systems because the non-linear cost function 15 
results in a 'direction of descent" (gradient) that is a functton of the parameters (variables) to be optimized. 
Neither projective transformation nor its variants (e.g., alfine scaling) rectify this situation. 

Summary of the Inventien 

Our Invention obtains a solution for the resource allocation task in a manner that, In part, is similar to that 20 
described by Vanderbel In U.S. Patent 4.744.026. Issued May 10. 1988. Specifically, the resource allocation task 
Is solved iteratlvely. and at each iteration a transformation is performed prior to affine scaling. The 
transformation rectifies the aforementioned problems by transforming the modeled system to convert the 
convex costfurtction (having ellipsoidal constant cost surfaces) to spherical form (havir^g spheroidal constant 
cost surfaces). Once transformed, affine scaDng is applied and th© iteration is continued. 25 

Because of this transformation, the search for an optimal point in the polytope becomes very efffclent due to 
special geometrical properties of spheres. Without this transformation, the search would be very inefRcient or 
in some cases it may even fail completely. 

Similar to what one may experience in the LP case, sometimes it is possible to get "attracted" near a wall of 
the polytope, which may not contain the optimal solution. This potential difficulty is overcome by "pushing" the 30 
"attracted* solutions away from the wrong wall, along the cunred cost surface, to a tjetter position. This 
procedure speeds up tr^e solution procedure considerably. 

As indicated above, the QP problem (quadrafic problem) Is a special case of the general NLFLC problem. 
V\/herea$ In the QP problem the cost function changes in a relatively slow and predicatable way. in the general 
HiFLC problem the characteristic of the cost surface can change dramatically from point to point. To obtain an 35 
improved cost allocation in the presence of such a varying cost function, a line search operation is performed 
to track the tecal optimal point from step to step to insure that a step Is not taken that passes beyond the 
minimum cost at the chosen direction. 
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FIG. 1 depicts the convexity property of a function; 

FIG. 2 is a graphical representation of the special case of optimal resource allocation for a QP problem ; 
fig! 3 Is a graphical representation of the general case of optimal resource allocation for a QP problem 
along with the iteration steps of the new algorithm; ^ 
FIG. 4 Is a depictk>n of the "zlg-zagglng' phenomenon of search directions due to ellipsoidal cost 

surfaces; 

FIG. 5 is an lllustratton of the computational advantage due to spherical cost surfaces; 
Rg! 6 is an illustration if the affine scaling concept; 

RG. 7 Is an illustratton of the potential push on curved cost surfaces for a QP case ; SO 
RG. 8 is a graphical representation of the concepts of the 'projected normal" vector and the "prolected 

push* vector used for the potential push; 
FIG. 9 is an illustration of the motBfication incorporated to the push used for QP, to cater the general 

NLPLjCcass; 

RG. 10 is a general flow chart which shows the procedure whteh separates the OP from general NLPLC 55 
problems and chooses the appropriate process to sohra them ; 

RG. 1 1 is a flow chart which describes the solution procedure ^process A) of a QP problem using the 
new algorithm; wi 

RG. 12 Is a flow chart wWch describes the solution procedure (process B) of a general IMLFLC problem 
(otherthan QP problem) using the new algorithm; and so 

RG. 13 Is a block diagram of a resource allocation/control system using the ideas described FIGS. 1 
through?. 
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DetaHed Description 
The Task 

The foOowfng description addresses the QP problem first becaus it simpler to understand and because it 
5 is an important subset of the NLFLC set. Thereafter, the results are extended to the g neral class of NLFLC 
problems, and some modifications are incorporated. Also, for the sake of simplicity and taking into account of 
the limitations of human perception, all of ttie graphical Blustratlons are confined to two dimensions only. 
IHowever. tha number of dimensions (number of decision variables) in practical problems can be hundreds, 
thousands, or even miUions, Moreover, without losing generality, only cost minimlzafion is discussed, but It Is 
10 understood that benefit maximization can be treated in a similar fashion. 

FIG. 2 presents a two dimensk>nal polytope 11 which encompasses the set of aD operafional states of the 
system. Every point within the polytope corresponds to a feasible state of allocation. Thus, point 14 Inside 
polytope 1 1 is a feasible point In FIG. 2. the elliptical contours 12 and 13 are constant cost surfaces. That is. 
(e.g.. state 15) the cost of allocation associated with aJi points on surface 13 is identical. 
IS Furthermore, the cost is smaller than the cost of allocation of all points on surface 12 (e.g., state 14). The 
challenge in connection with the FIG. 2 system is to arrive at a state which results in the least allocation cost 
(conresponding to the inner most ellipse) and at the same time maintain feasibility (stay inside the polytope). 
This arrive-at state is the optimal solution to the allocation task, which In FIG. 2 is point 16. It may be noted that 
point 16 is the optimal allocation with no constraints considered. This is a very special case of the QP problem 
20 and such occurrence is extremely rare In real-life physical resource allocation problems. 

A more realistic QP case is illustrated by FIG. 3. The unconstrained optimum 21 is not feasible since it is 
outside polytope 24. Therefore, the optimal state for the FIG. 3 system is point 22. which lies on wall 23 of 
polytope 24. The optimal cost Is the cost associated with the ellipse 25. 
Our Invention starts with an initial feasible point 26, and proceeds radially, but always staying in the interior of 
2S polytope 24 (thus maintaining the feasibility), in steps 29, 30. etc.. to states 27, 28. etc., each closer to the 
optimum feasible state 22. 

The QP Slttjatfori 

In standard vector notations, a typical QP problem can be expressed as follows: 
30 minimize 

J- Vax^Qx-c^x (1) 

subject to a set of linear constraints 

Ax » b (21 

and a positivity constraint x ^ 0 on x where 
3S X is an n component vector called the "solution vector* 

Q is an n X n matrix (at least positive semi-definite) called the *cost matrix' 
c is an n component vector called "cost vector" 
A is an m X n matrbc called the "constraint matrix* 
and 

40 b is an m component vector which depicts the constraint limits. 

The above representation of the QP problem is commonly refenned to as the "standard fonm*. See Gill et al.. 
referenced infra. All objective functions and ail constraint relattonships pertaining to QP proglems can be 
reduced to this form by simple algebraic manipulattons, with the help of what are known .as 'stack* and 
"surplus" variables. These technk)ues are well-known in the prior art, as described, for example, in the 
45 afbrementtoned 4,744.025 Vanderbel patent. 

For the special QP case depicted In FIG. 2 (where the alasolute minimum is within the polytope) tha optimal 
solution and the corresponding cost can be easily catoulated as, 

Xoirtlmian = OT'^C (3) 

and it is straightfonivard to compute 

SO 



(4) 

55 

For the general case of QP problems (where the unconstraint optimum lies outside the polytope) the 
solution can proceed by following the gradient of the cost very much as it is done by Karmarkar. However, this 
procedure may be Inefficient for the reasons explained below. 
With reference to FIG. 4. lines 32, and 33 are the gradients (hence the search directions) which are 

60 orthogonal to constant cost surfaces 35. and 36. Proceeding in accordance with gradient concept alone, the 
solution would start from one state, say 30, and proceed along th gradient (32) to another state, say 31 . State 
31 is selected based on the criterion that Is related to the change in the gradi nt; and more specificalty. based 
on a determination that the current direction (32) Is tangential to another constant cost surface, such as 
surface 36. When that condition is detected, a new gradient Is computed and the solutton proceeds from state 

6S 31 to the next state, such as th state represented by point 48. The aim is to reach the optimal state 34 in this 
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manner. But. as can be observed, the procedure "zig-zags' considerably before reaching tWs point. 

Our invention overcomes this difficulty by transforming the elliptical (or eilipsoidaMn multi-dimensions) 
cost contours to circular (or spheroidal- In multi-dimensions) contours. Thus, with reference to FIG. 5. 
because of the special geometrical properties of sph roidaJ surfaces, the search direction 41 Is always 
orthogonal to contours 42. 43, 44 etc.. and this search dir ctton is oriented closely toward the optimal state. 5 
45. 

The mathematical translation of the above idea Is the transformation of the variable x (solution vector) to 
x' - L'x (5) 

where L is a vector that Is related to the cost malrtx Q fciy 
Q^UJ. (6) 

Note that X ^ 0. but x* is unconstrained. 

Following this transfornr^tion. the entire problem is scaled (affine scaling) by mapping the cunrent state 
(point) to a common distance ji from the constraint walls, where |i is a poslthre scalar. This process can be 
mathematicaUy represented by the product ir'x. where 5 is a diagonal matrix YffXU diagonal elements set to 
(1/^)xi- Refening to FIG. 6 where the unsealed and scaled polytopes are shown, the current stateSS In the is 
original (unsealed) poiytope 51 is mapped to state 56 in the scaled polytope 52, forcing it to be at a distance (l 
from walls 54 In the positive quadrant 

With some mathematical manipulations, it can be shown that the negative of the direction of movement 
(step direction) from a state to minimize the cost of allocation Is expressed by. 

5iCp=--{i-HAT(AMAn-^AlH(QxO-C) (7) 20 
where 

H equals (XQ + IV»'*>"^. 

Dxo is a diagonal matrix with the diagonal entries the components of x<>. and 
X. is a positive 'scale' parameter, which advantageously is set to ir^. 

The next task Is to calculate the step length which determines the size of the step (movement) away from x» 2S 
in the direction 5xp. 

One constraint, of course, is the presence of the walls of the polytope. The step length due to this constraint 
Is limited by 



„,. (8) 



max. 



[8«pT(Q»°-c)] 



oa = : .... _ . (9) 



30 



which goes takes up to 97QAd percent of the distance to the nearest constraint watl. Another constraint. 35 
however, relates to the curvature of the cost surface. As shown in FIG. 4. while proceeding In the direction Sxp. 
a point is reached where the cost stops to decrease and tiegins to increase. Going beyond this point is<:lear1y 
not advisable, and therefore, another value is computed that evaluates the step length to the minimum cost 
along the chosen cfirection. This step length is 



40 



45 



The actual step length ts, of course, the smaDer of the two. i.e., 

a « mln(ai,a2) . . 

The new state x^ is obtained by evaluating. 

xi-xa-o6xp. (10) ^ 

The expression for the descent direction. 8xp. is called the 'null space" fonnulation. because it is obtained by so 

the proiectton of the gradient {Cafi - c) on to the null space of A. An alternate fonn of this would be the 'range 

space* formulatkxi. The only difference between the two formulations is the way in which the projection matrix 

Is constmcted. it can be shown that the direction for the range space fbrnnulation of the QP algorithm Is given 

by 

6xp«-r(ZXQ + irM2n-'Z(Cbt©-c). (11) ^ ^ 

where Z Is cafled the 'orthogonal compUmenf matrix of A and satisfies the feasibility condition ATv = 0 for 
any vector v. 

To construct Z. the matrix A may be represented by [ B I N ] where B and N are of dimenstons m x m and 
mx(n-m) respectively, and then 

Z = [(-B-^M)» I O (12) ^ 
Where I represents the "identity matrix" of dimension m. The scale parameter X can be adjusted to make the 
conv rgence of the algorithm faster. (It can be shown that for the LP case, the scale is totally absorbed In the 
step length ai , and consequently scale adjustments are Ineffective). The way to do this is to perform a 'search 
on X, but it should be noted that every time th value of X Is changed, a matrix inversion must be performed. 
Therefore, the following recursion is useful In adjusting X in the nelght>orhood of some X** without performing 65 
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the inversion each and every time. 

Sxp^'xa'-xy-g^*) (13) 

5 s=0 
where 

fO po Is the projection matrix. [I - HA^CAHAVAJH. evaluated at X^, and 

9« is the gradient of the cost function evaluated at ifi which in this case is Cbc® -c. 
Also, sometimes it is possible that the successive solutions of the algorithm to be near a constraint wall. 

Under that situation, the convergence of the algorithm can be slow. When this happens, a potenttel search 

method can be used to recenter the solution on the current cost surface. With reference to FIG. 7, 61 is the 
fS polytope and 62 is the state from which a move is made in the direction 63 to reach the state 64, which is near 

wall 65. The optimum is state 66. What is necessary to do is to slide 64 along the constant cost traiectory 67. to 

recenter at a better position, such as position 68. This Is what we call the potential search. 
The essential idea beNnd the potential search Is to construct a 'push vector' which pushes the solution 

away from the walls. Of course, the path traced by this push vector should be of constant cost, while 
20 maintaining feaslbllit/. 

The push Is accompUshed by evaluating tip which is the projected normal at xo. and by evaluating Vp which is 

the projected push direction at X0, where xo Is the point that is to be pushed away from the walls. Letting « be 

the path traced by xo so that the cost is unchanged, we compute 

Hp = ^(ZQZVZgo (14) 
^ and 

Vp - C r(ZHr)-i2vo (15) 

where 

(i) Z Is such that A2^r - 0 for any vector v. 

(ii) go Is the gradient to the cunred level set at x^, 

30 m) 



35 




Vo'^Z^(ZHZ'''rlZQZ^(ZHZ'^r^ZvoJ ' 



(lv)H = (XQ -hD-^.and 
(v) vo Is the "Karmarkar Push Vector' obtained as 
40 voi »xoi''^-Pgoi fori « 1.2, --n 

where B = — — — and 
gogp 

^ v'p-C2^(ZH2f)-^Zxo-\gp-C2T{ZHr)-^2go. 

Letthg a be HQIi, where || B denotes the spectral norm, a "search" along the trajectory 
xt » xo + (cos at-1)np + (sin a t)Vp (16) 
can be performed so that the 'potential function" 

-Flog(x^) 

SS Is minimized. Referring to FIG. 8 as an example. 71 depicts the spherical surface. 72 and 73 are the projected 
nonnal and the projected push vectors at state 75 vectors, and curve 74 represems the trajectory of constant 
cost. 

Generalizing to NLFLC 

so Our QP algorithm can be generalized to other NLFLC optimization problems In systems where the cost 
function is other than quadratic but stOI convex. The task Is then represented by 
min f(x) (17) 

subject to a set of linear constraints 
Ax » b and x2sO. (18) 

65 Letting g(x) be the gradient of f(x) at x. (i.e., VM = flW). and Gfx) be the Hessian matrix of f at x (l.e., 
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72f(K) = G(x)). and replacing Q by G(x) and Ox - c with gM In our QP algorithm, yi Ids a new atgornhm to 
solvB the NLFLC problem. Unlike the QP algorithm, however, the Q term depends on the vector x, and this 
InvaJldates our analytical xpresslon (9) for the step size 02. The step size 02 In tNs case Is then determined bi 
a fast 'line search method', where th step length of desc nt Is determined by successive evaluations of f(x) 
in th descent direction. ^ 

incorporating the required ilne search, otir new algorithm for solving the NLFLC problem can be stated as 
follows: 

Choose two positive constants p and y such that 0 < P < 1 and 0 < y < 1. along with a starting feasible 
solution Then, compute the search direction 6xp from equation (7) or (1 1 ) by substituting G(x«0 and g{x^) In 
place of Q and Qx^} - c respectively. Use 02 In (9) as an estimate for the potential step size, where Q is replaced 10 
by G(xO). Lastly, obtain the true step length a such that 
f(x«-a€xp) ^ l(xO)-YagT(x^Xp (19> 
where 

a - m^n(a2P^a1) (20) 

the term 01 is given by equation (8) and 6 is the first positive Integer (Including zero) that satisfies equation IS 
(19), it can be shown that the Ilne search algorithm converges (reaches the optimal answer) under reasonably 
weak assumptions. 

The rest of the procedure Is the same as the QP algorithm, except for the following difference in the 
potential search. While perfonning the potential search, we recenter along a path based on a quadratic 
approximation of f(x). This recentering scheme maintains feasibility but does not keep f(x) constant. 20 

Refening to FIG. 9 for example, state 81 is translated In direction 82 to state 83. This solution is then pushed 
along the constant cost surface 85 based on the quadratic approximation of f(x). The actual constant cost 
surface is 84. and we stop the push along the trajectory of the approximation of f(x) when the difference 
between trajectories 84 and 85 (which In this case Is distance 86) exceeds a predefined value, say 5J. 
TTiereupon, we start a new iteration. ^ 

Finally, a stopping criterion for the new algorithm (for QP as well as NL5LC) could be 
8xp S e (21) 
where e Is some fbced small positive number and 

p.A^(AHAV AH] [QX9-C] g 0 (22) 30 

The Process ^ * 

The prooess of Improving the operational state of a commercial enteipiise such as Improving the state of 
resource allocations of the enterprise Invokes first Identifying the resources that are available, determining the 
resource allocation needs and constraints, and selecting a cost function to be minimized that makes sense for 
the particular circumstance of the enterprise. All of this Information Is then represented by symbols (variables) 35 
for easy manlpillation. and the task is stmclured in terms of the represented variables as a minimization task. 
When so structured and when it Is found that the constraint relationships are linear but the cost function is 
convex, then the prindpies of our inventton can be applied to an advantage. 

In such a circumstance, according to FIG. 10 which outlines our method in broad terms, the task Is 
formulated as a QP or general NLFLC task in block 91, and In block 92 an initial feasible state of the system is 40 
selected <x<^, and a scale parameter \ is chosen (usually, 1.0). This Is called the "phase 1' problem. As the 
constraints are linear, one can use the same phase 1 operatton of LP case. This process is dlsctosed in the 
aforementfoned U.S. Patent No. 4.744.02B issued to R. J. Vanderbel. Of course. In existing systems the task of 
chosing an Initial state is trivial since the existing state of the enterprise can be selected as the Initial state. 
Continuing with RG. 10. btock 93 represents a process to detect whether the optimization problem Is QP or 45 
general NLFLC. H It Is QP. process A (btock 94) is Invoked; othewvise process B (block 95) Is Invoked. Thus, 
block 96 represents the entire iterative scheme (though In very general tenns) to solve all NLFLC problems. 
Lastly step 97 makes the assignment, as determined by block 96. to place the enterprise in the neighbortiood 
of the optimum stale of the system. Processes A and B are described In more detail In FIGS. 11 and 12. 

Referring to FIG. 1 1 . block 101 is the detection step for determining whether the number of variables (n) Is SO 
more than the number of constraints (m). If (m > n). then we proceed to btock 103 to use the null space 
formulation to obtain the descent direction. OthenMse. we go to block 102 to use the range space formulation. 
Next, the two step lengths (at and 02) are computed In btock 104. and block 105 performs the translation from 
yfi to in the direction of descent using the step length computed above. After that, block 106 determines 
whether a scale adjustment of the initially chosen scale vahie should be done, and If advisabte, the process S5 
continues to block 1 07 where the adjustment is made, a new direction is evakrated. and control is returned to 
block 105. When it is determined that no scale adjustment Is necessary. control transfers to block 108 where. If 
necessary (that is, if a = at) . a potential push Is perfomied in block 108 to recenter the solution, as described 
before. Finally, block 109 checks the convergence condition, such as the one expressed by equations (21) and 
(22). When the condition is met. process A terminates; otherartse control returns to block 101. The €0 
above-described steps are repeated until the convergence oondltlons are met. 

Process B is very similar to process A and it is illustrated in FIG. ia Block 110 perfonms a quadratic 
approximation on f(x). Blocks 111, 112, and 113 are the same as blocks 101. 103. and 102 or process A. In 
block 114 we compute th st p I ngth as defined by equatton (19) and use it Inconjunctton with th search 
direction ttiat Is obtained from block 112 or 113. This op ration is done In block 115, whteh essentially uses & 
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equation (20). Blocks 116 and 117 are Identical to blocks 106 and 107 of process A. In block 118. we pertorm 
the approxiniat potential push which was described b f re in conn ctfon with FIG. 9. essentially using 
equations (14) through (1 6). Block 1 19 is identical to block 109 and the whole process B Is repeated until these 
conditions are satisfied. 

5 

The Apparatus 

FIG. 13 illustrates a process control system in accordance with the principles of our invention. It nrtay be a 
telephone communication system, a manufacturing process, a navigation process, a financial system (e.g., 
portfolio management system), a lease of planes (as described before) arrangement, or any other industrial. 
10 technotoglcal or economic process which is to be optimteed and which exhibits the 'convexity' property 
discussed above. 

A cost function evaluation device 170 receives the cost data (information on how allocation cost varies) on 
bus 171 from a computer terminal or some other separate process. Device 170 aims to accommodate the fact 
that costs are not necessarily fixed and that different enterprises may choose different cost arrangements that 

IS they may wish to minimize. Device 170 fits the incoming information to a cost function representation, which 
for the systems under consideration is in general non-tlnear and more particularly, convex. The information is 
passed by device 170 to cost evaluator 172 which evaluates the value of the fitted cost functton at any state of 
the system. That is. given an operational state of the system represented, device 172 evaluate the cost In 
accordance with the given cost function. To evaluate the cost, device 172 must be provided with a description 

20 of the state of the system, and that information is provided by controller 185 via lead 173. The cost 
determination is passed to controller 180 via lead 174. 

Similarly, a limit register 183 is provided to store the physical limits on the allocation variables. These Omits 
can be eritered into the system via lead 184. like the cost data. 
The system of FIG 13 also Includes sensors 186, 187, • • 188 whk:h dynamically sense the constraint 

25 coefficients of the enterprise. In the airiine industry, as an example, the number of planes at different locations 
win. of course, change with time, and even the number of operating airports will fluctuate. It Is tfiese 
'environmental" changes that sensors 186-1B8 monitor. Each constraint sensor 186-188 has a correspondirrg 
change detector, 189. 190, • • • . 191. The detectors detect the changes in the output of each of the respective 
sensors. A change indication signal from each of these detectors is applied to change bus 192 and thence to 

30 AND gate 193. Also applied to gate 193, on lead 194, is a signal from controller 185 indicating the termination of 
the optimization procedure. Through a separate oul^^ut port at each of the change detectors, the outputs from 
sensors 186-188 are applied to controller 185. 

Controlter 185 of FIG. 13 Is. In the preferred embodiment, a programmed digital computer having stored 
therein the program which implements the ffewrcharts of FIGS. 10. 11 and 12, making use of the ideas 

35 presented in RGS. 1 through 9. Controller 185 may also comprise a complex of tiardwlred circuits, designed to 
cany out the operations described by the flowcharts of FIGS. 10 throu^ 12, or to realize the kieas of FKaS. 1 
through 8. 

Since NLFLC controller 165 of FIG. 12 utilizes the extremely fast procedures illustrated in FIGS. 3. 4. or 5, 
control values are availatsle for registers 195-197 in a very short time. Moreover, as the constraints change. 
40 these changes are sensed by sensors 186-188, detected detectors 189-191, and used partially to enable 
AND gate, 193. When the procedures of FIGS. 10-12 are complete, controller 185 generates process control 
signals and transfers the signals to registers 195-197. Simultaneously, controller 125 generates an enabling 
signal on lead 194 to AND gate 193, completing the enablement of AND gate 193. The entire process is then 
repeated. 

45 A typical type of problem in modem control theory area to which the present invention can be applied is 
stochastic control of systems, for example the control of robots. A description of these type of problems can 
be found in the book "Computer Controlled Systems' by K. J. Astrom and B. Wittenmark, published by 
Prentice Hall In 1986. 

Other problems which would benefit from our invention are financial portfolio management, optimal loading 
50 problems of power utility companies, and optimal staff allocation problems. 

It should be noted that the matrices involved in most practical NLFLC cases are sparse in nature and sparse 
matrix techniques can be used in the operations depicted In RGS. 10 through 12. 



55 Claims 

1 . A method for allocating industrial facilities in an enterprise so as to reduce the cost associated with a 
selected allocation of said facilities In said enterprise, wtrere said allocation of said facilities is linearly 
constrained by a set of constraints and said cost is expressed in terms of a given convex function, said 

60 method comprising facilities allocation iterations wtiere each Iteration receives a potential facilities 

allocation having an associated cost,and develops a potential facilities allocation for ttie next iteration as 
long as a selected stopping criterion Is not met. each iteration further comprising the steps of: 
performing a scaled spheroidal transformation that op rates on representation of said received potential 
facilities allocation to form a transformed repr s ntatlon of said potential faculties allocation such that 

65 . components of said transfonned representation of said potential facilities allocation ar equal laffin 
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scaling), and th sum of squared components of transformed representations of other facilities 
allocations having the same associated cost is equal to the sum of squared components of said potential 
facflitiss allocation (spheroidal transfonnation}: 

developing a translation vector for use In creating said new potential faculties allocation for the next 

tes«^he components of said translation vector to develop a process-stop signal when said translation 
vector fans below a preselected ttveshotd: 

developing said new potential facilities aDocation for the next iteration when said translation vector fails to 
fall below said preselected threshold: and 

applying the last developed new potenttal facMitles allocation to said enterprise when said translation 10 
vector faDs below a preselected threshold. 

2. A method for allocating Industrtal facilities In an enterprise so as to reduce the cost associated with a 
selected allocation of said facttities. where said allocation of said facilities is linearly constrained and said 
cost Is expressed In terms of a given convex function, said method comprising facilities allocation 
Iterations where each Iteration receives a potential facilities allocation, and develops a potential faciKties 15 
allocation for the next iteration as long as a selected stopping criterion is not met, each Iteration further 
comprising the steps of: 

transforming the representation of the potential facilities allocation, vector to develop a spheroidal 
cost function for said enterprise such that the sum of squared components of the transformed facilities 
allocation representation (components of x^) is a constant value and the sum of squared components of 20 
representations of transformed other faciiitalos allocations (components of x^'^x^, - • 0 having the same 
cost as the cost of said potential transfonmed facilities allocation (x*) is equal to said constant value; 
scaling the representation of the potential facilities allocation, x». to move the potential facilities allocation 
to a seated potential facilities aUocation represented by vector x«" and having identical scaled facility 
allocations (components of vector x"*) ; ^ 
developing a translation vector for use In creating said new potential facilities allocation for the next 

iteration; . , u i 4- 

testing the components of said translation vector to develop a process-stop signal when said translation 

vector falls below a preselected threshold; 

developing said new potential facilities allocation for the next iteration when said translation vector fails to 30 
fail below said preselected threshold ; and . ^ . . 

applying the last developed new potential facilities allocation to said enterprise when said translation 
vector falls below a preselected threshold . 

3. TTie method of claim 1 where said step of performing a scaled spheroidal transformation comprises 

the steps of 

transforming the representation of the potential facifitles allocation to develop a spheroidal cost function 
such that the sum of squared components of the transfonmed facilities allocation representation is a 
constant value and the sum of squared components of representations of transformed other fadtilies 
allocations having the same cost as the cost of said potential transformed facilities allocation is equal to 
said constant value, followed by the step of ^ 
scaling the representation of the potential facilities allocation. x*». to move the potential facilities allocation 
to a scaled potential facilities allocation and having identical scaled facility allocations. 

4. The method of claim 1 where said step of performing a scaled spheroidal transformation compnses 
ttie steos of 

scalfrig the representation of the potential facilities allocation. x». to move the potential facilities allocation 4S 
to a scaled potential facilities allocation and having Identical scaled facPity allocations, followed by the 

toOTsformlng the representation of the potential facilities allocation to develop a spheroidal cost function 
such that the sum of squared components of the transfomned facifities aUocatlon representation is a 
constant value and the sum of squared components of representations of tnansfomned other facilities SO 
allocations having the same cost as the cost of said potential transformed facilities aUocation is equal to. 

said cx)nstant value. ^ * ^ ^# 

5. A resource aUocation system for allocating facilities of an enterprise so as to reduce the cost of 
operating said enterprise with a selected aflocation of said facilities, where said allocation of said facilities 

Is llnearty constrained and said cost is expressed In terms of a given convex function, compnsing : 55 
an N plurality of sensors, where N Is a positive integer, for determining the existing state of allocations of 
safdfacilities; 

means for applying an allocation of said fadntles; 

means responsive to said sensore for storing output signals of said sensors in an alterable memoiy, 
thereby forming a stored state of allocallons of said facilities, said state having N components; «? 
transform means responsh^e to said means for storing for performing a scaled spheroidal transfomiation 
of said stored state of allocations to form a transfonned aiocatlon state signal such that components of 
said transformed allocation state signal are equal, and the sum of squared components of transfonned 
repr sentatlons of other facilltales allocation states having the same assodat d cost is equal to the sum 
of squared components of said transf omied allocation state signal ; 
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translation menas responsive to said transform means for developing a translation vector signal; 
comparison means responsive to said translation means for developing a process-stop signal when said 
translation vector signal falls below a prsselected threshold ; 

means responsive to said comparison means for developing an improved fadfittes allocation state and 
storing said improved facilities allocation state in said alterable memory: and 

means responsive to said process-stop signal for applying most recent one of said Improved facilities 
allocation states to said means for applying an allocation of said faciiities. 

6. A system for optimizing the performance of oontroBed process in accordance with an optimization 
criterion, said system comprising: 

process control devices for controOIng said process In response to a control signal set through control 
signals; 

a plurafity of sensors for sensng variable conditions affecting the operation of said process via input 
signals; 

a non-linear programming controller interposed between said sensors bikJ said process control devices 
for developing control signals applied to said process control devices, where the developing of said 
control signals is done in accordance with said optimization criterion and where said optimization 
criterion is a non-linear convex function of said input signals ; 

said controller including means for iteratlvely Identifying successive tentative strictly feasible control 
signal sets and selecting each tentative control signal set in the directton of the steepest gradient of a 
normalized version of the optimization criterion in a space where the constraints on the control signals are 
alsonomnalized. 

7. A method for improving the operational state of a commercial enterprise characterized by a number 
of resounces D. (i = 1 .n) among a plurality of physical resource users subject to constraints Aipu » b| and 
xi^ 0(1 1.n;j » 1.m) in such manner to optimize a cost function (l/2)XiQQXi-cixi(i » 1.n),said 
method comprising the steps of: 

(a) selecting an Initial allocation ifi » - - ^ meeting the said constraints. 

(b) using the transformation « L^x where 

Q B LL^ to make the constant cost surfaces spherical, 

(c) using the scaling y » D'^x' wttere 

D » dlag(1/^Kx?>4 ' - to place x^ at a distance (x from the constraint walls. 

(d) computing the descent direction 5xp by using the relationship 
5xp » -[I - HA^(AHAV A]H(QxO * c) where 

H » (Xa + Dxo-2)-i,Dx9 is the diagonal matrix containing the component8orx<^» 

X = 11-2, and 

I is the identity matrix, or 

6xp = -Z^CZCXQ + D-Mi^-*Z{QxO - C). 

where Z = ((-B-^Nr 0. where matrix A <- [B M] 

(e) computing the step length a mln (ai>a2), where 



(f ) selecting the new allocation x^ by x^ ^ x^ - a5Xp, 

(g) adjusting X for maximum cost reduction, 

(h) recentering In the new constant cost surface If a « ai using the potential push described in 
the detailed descr^tion of the algorithm, 

(i) testing for the criterion 5xpi g e (where 6 Is iwme fixed small positive number). 
(J) testing the criterion [1 - A^CAHAV^ AH] [QX<> - cj ^ 0 

(k) stopping the iterative scheme if (I) and (]) are satisfied, or 

(I) going back to step (d). replacing x^ by x^ are repeating, 

(m) allocating said resources In accordance with x^ If (i) is satisfied. 
8. A method for allocating industrial, financial, or technological resources xi (i » 1.n) aniong a plurality 
of physical resource users subject to constraints Aqxi » b| arxl xi ^ 0 (i = 1 .n,j « 1 .m) in such manner to 
optimize a general convex non-linear cost function f(x), said method comprisintg the steps of: 

(a) making a quadratic approximation on f(x) and performing the steps of (a) through (e) of claim 
10. 

(b) calculating the step length a as per the lin search rule 
Hyfi - a5xp) ^ f(xO) - yag^(xO)6xp wh re 

a « min (aaP^.ai) where h is the firsl positiv integ r (including zero) that satisfies the said line 
search rule, 

(c| obtaining the new allocation as per (f) of dalm 10, 
(d) periormingascaleadjusmentasin (g)ofGlaim 10, 
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( ) recentering on constant cost surface as of (h) of claim 10. with the modification that the 
translation should k>e stopped, also taking into account the difference tietween the actual functional 
value and its quadratic approximation. 

(g) allocating the resources Itlce (m) of daim 10. 

9. Apparatus for allocating resources in an optimal manner among users of said resources comprising 5 
means for receiving information concerning said users, cone ming the availability of said resources, and 
conceming the constraints on the allocation of said resources. 

means. utOirlng the algorithm of claffn 7 or daim 8. for iteratively approximating the optimal allocation of 
said resources among said users. 

means for allocating said resources in acconjance with the last option of allocating said resources 10 
approximated by said iterative approximating means. 

10. The method of claim 1 where said step of developing said new potential facilities allocation for the 
next iteration comprises a line search along the direction of said translation vector. 

11. The method of claim 1 where said step of developing said new potential facilities allocation for the 

next iteration comprises a line search along the direction of said translation vector to identify a potential l$ 
facilities allocation whose cost is not higher than the -cost of a potential toilities allocation further along 
said translation vector direction. 

12. The method of claim 1 where said step of developing said new potential facilities allocation for the 
next iteration comprises a Dne search along the direction of seud translation vector to identify a potential 
facilities allocation whose cost is not higher than the cost of a potential facilities allocation furttter along 20 
said translation vector direction and which Is not closer than a preselected threshold to violating any one 

of said constraints. 
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